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Abstract 

We study a simple replicator model with non-symmetric and nonlocal inter- 
actions. Hierarchical structures with prey-predator relations are self-organized 
from a homogeneous state, induced by the dynamical instability of nonlinear 
interactions. 
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Ecological systems and evolution dynamics have been studied as examples 
of complex systems PE]. The replicator models are often used to study bio- 
logical evolution among mutually interacting species [3]. The time evolution of 
population x\ of the ith species obeys an equation 

^ = %i(fiO + X! w hJ x J ~ c o)> for i = 1, • • ■ , JV (1) 

3 

where N is the total number of species, /io is the natural growth rate of the 
ith species, Wij denotes the strength of the interaction from the jth species 
to the ith species, /$ = fio + Y) j w i,j x j is the growth rate for the ith species, 

and Co = Yli=i fi x i/J2iLi x i IS the average value of the growth rate. The sum 
of Xi is conserved in the time evolution of Eq. (1). If the sum is assumed to 
be 1, Xi has a meaning of the population ratio of the ith species. The above 
equation implies that the populations of the species which have the growth rate 
larger than the average value increase. If Wij 's are all zero, only the species 
with the largest natural growth rate survives. It implies the survival of the 
fittest species. If Wij's are not all zero, the coexistence of several species is 
possible. Complex dynamics and statistical properties of the replicator models 
with random interactions have been studied 0]. Steady evolution called "Red 
Queen" evolution can be modeled with this type of equation 0- We consider 
the case of fio =const. in this work, that is, all species are equivalent by nature 
in the evolution dynamics. 

The population dynamics using Eq. (1) is closely related with the Lotka- 
Volterra type equations, which have been more intensively studied in ecological 
systems 6 . In the usual Lotka-Volterra models, the interaction strength {wjjjj-'s 
are usually assumed to be constant. They are usually non-symmetric, that is, 
w i,j 7^ w j,i- ln the prey-predator relation, where the jth species is a prey of 
the ith species and the ith species is a predator of the jth species, Wij > 
and Wj t i < 0. In the competitive relation between two species with equivalent 
strength, the interactions may be symmetric as Wij = Wjj < 0. However, we 
can consider some cases where the interaction strengths change dynamically. For 
example, there are some models, in which Wi j is a function of the populations 
Xi and Xj. In the population dynamics of trees in forests, competition between 
trees is not symmetric in general. Larger trees prevent the growth of lower trees, 
but the lower trees do not seriously affect the larger trees. Assuming that x.- L is 
the height of a tree at the ith position, Yokozawa and Hara proposed a model, 
in which W(xi,Xj) = expressed as W(r) — r 13 for < r < R c and 

W(r) — R@ for r > R c , where r = Xj/xi 0. We proposed a similar form of 
nonlinear interactions expressed as 

Wij = 0, for r < di, 

= ai(r — di), for d\ < r < d2, 

= a 2 (r - d Q ), for d 2 < r < d 3 , 

= a^(r — di), for da < r < d^, 

= 0, for d 4 > r, (2) 
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where r — Xj/xi, d\ = l/t&t, d 2 = l/cfa, d\ < d 2 < d < 1 < d 3 < dj, 
and ai > 0, a2 < 0, 03 > 8 . The interaction Wij(r) is a continuous and 
piecewise linear function of the ratio r = Xj/xi- This form of interaction seems 
to be peculiar, but it is one of the simplest forms of continuous functions, by 
which both the mutually equivalent competitive relation and the prey-predator 
relation can be described. The interaction is competitive and equivalent, i.e., 
Wij = Wj t i = 02(1 — do) < 0, when the population sizes of the ith and jth 
species are equal. Even if the population sizes between the two species are not 
completely equal but are nearly equal {do < r < (1/do)), the interactions are 
competitive, i.e., Wij < 0, wj^ < 0. If the ratio r = Xj/xi of the population 
sizes is even smaller and satisfies d\ < r < d , the jth species becomes a prey 
of the ith species (that is, the ith species is a predator of the jth species), i.e., 
Wij > 0, Wj t i < 0. If the ratio r of the population sizes is even larger and satisfies 
(1/do) < t < di, the ith species becomes a prey of the jth species (that is, the 
jth species is a predator of the ith species), i.e., Wij < 0, Wj t i > 0. If the ratio 
r is too large (r > d.4) or too small (r < d\), the interaction becomes negligible, 
i.e., Wi_j = Wj t i = 0. The ecological relations change dynamically with the ratio 
of the population sizes in our model ecological system. Parameter values of 
d = 0.85, dt = 3/5, d 2 = 2/3, d 3 = 3/2, rf 4 = 5/3, a 2 = —1.3, a x = a 2 (d 2 — 
do)/{d 2 — g?i) andd3 = a^dQ—d^/^d^—d^) are used in the numerical simulations. 
A homogeneous solution x< = 1/N for every i is a solution of Eq. (1), owing 
to the symmetry of the coupled equations. However, the homogeneous solution 
is not always stable. The present author, and Yokozawa and Hara found that 
layered structures are self-organized in a globally coupled model and a locally 
coupled model with the nearest and next-nearest neighbor interactions, as a 
result of the instability of the homogeneous solution |Hj . 

In this paper, we study a nonlocally coupled model based on Eq. (1) with 
nonlinear interactions of the form Eq. (2). The position of the ith element is 
located in a one-dimensional lattice or randomly distributed in two dimensions. 
We assume further that the interaction Wjj takes the value of Eq. (2) only if the 
distance between the ith and jth elements is smaller than a value determined 
by a function f(x), where x is the average population size {xi + Xj)/2. We 
use a simple function fix) = ax in numerical simulations. The function f(x) 
expresses the range of the interaction, and the function f(x) — ax implies that 
the interaction range increases in proportion to the population size. 

Firstly, we study a one-dimensional model. Each species is located at i, 
where 1 < i < N. The homogeneous solution = l/N. The interaction 

strength Wij takes the value of 02(1 — do) for a pair of the same population size. 
The linearized equation for a small perturbation 8xi = X4 — x s around the 
homogeneous solution obeys 

dSx 

= x s{^2 ° 2 ( 2 _ do)8xj - a 2 <fai}) for i = 1, • • • , AT, (3) 
j 3 

where the sum is taken for the range satisfying |j — i\ < a/N . The eigenvalue of 
the linearized equation can take positive values for a/N > 1. For example, if the 
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site number satisfying \j — i\ < a/N is 5, the eigenvalue Afc for the wavenumbcr 
k is expressed as = x s a2{(2 — dp)(l + 2cos/c + 2cos2fc) — 5}, and therefore 
A fe > for cos -1 { (-1+^5 + 20/(2 - d ))/4} < k <ir. Under the condition, the 
homogeneous state becomes unstable. We have performed numerical simulations 
for N — 200. The initial condition was Xi — \/N(\ + 5xi(Q)), where Sxi(0) was 
chosen as a random number between -0.01 and 0.01. Figure 1 displays the 
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Figure 1: Profiles of population sizes x% for a one-dimensional model Eq. (1) 
with N = 200 and (a) a = 180, (b) a = 300, (c) a = 1000, (d)a = 4000. Solid 
lines denote pairs with prey-predator relations. 



population size x 4 for (a) a = 180, (b) a = 300, (c) a = 1000, (d) a = 4000 
at t = 20000. Solid lines denote pairs with prey-predator relations, in which 
Wij < and Wj : i > are satisfied. The homogeneous state is stable for a = 180, 
since there are no interactions except for the self interaction expressed by Wi^. 
The homogeneous state becomes unstable for a/N > 1 and zigzag structures 
appear for a — 300. As a is increased, hierarchical structures with higher 
levels are created. Strong elements appear and dominate locally the regions 
of their population sizes. The local winners have larger population sizes and 
therefore larger interaction ranges. Competitive relations are therefore created 
among the local winners, and the higher hierarchical relations are created by 
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the competition among the local winners like a tournament. For a = 4000, 
there appears an element which dominates the whole region, and a hierarchical 
structure from top to bottom is created. Which element dominates the whole 
region depends on the initial conditions. 

Next, we consider a two-dimensional model. The position of the ith element 
is randomly distributed inside a circle of radius 10. The population size at the 
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Figure 2: Hierarchical structure for Eq. (1) with TV = 400 in two dimensions. 
Each element is randomly distributed in a circle of radius 10. The parameter 
value of a = 800. Figure 2(a) is the network with larger size populations, which 
satisfies Xi > 0.005, and (b) is the network with the smaller size populations. 
The area of the square at each site is proportional to X{. Figure 2(c) is a 
histogram of the population size. 

ith element is denoted as Xi. The system size is N = 400. The initial conditions 
are x$ = 1/N(1 + Sxi(0)), where 6xi(0) is a random number between -0.01 and 
0.01. As a result of the time evolution by Eq. (1), a two-dimensional hierarchical 
structure is naturally created. As the parameter value a is increased, we have 
observed development of hierarchical structures as in the one-dimensional model. 
Figure 2 displays a developed hierarchical network for a is 800. The prey- 
predator relations are displayed in this figure. To show the network clearly, we 
have classified the set of i into two groups. Figure 2(a) is a network for larger size 
populations, which satisfies Xi > 0.005, and Fig. 2(b) is a network for the smaller 
size populations, which satisfies Xi < 0.005. The population size Xi is expressed 
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by the area of the square located at the ith position. The elements, which have 
larger populations, dominate the larger regions. The average number of preys 
for each element is 4.3 in this simulation. Figure 2(c) displays a histogram of 
population size. It is seen from this histogram that a layered structure with 7 
levels is created. 

To summarize, we have proposed a nonlocally coupled replicator model for 
the self-organization of hierarchical structures. The interaction range was as- 
sumed to be proportional to the population size, and the nonlinear interactions 
between two species are assumed to depend on the ratio of the population 
sizes. Owing to the two assumptions, a hierarchical structure of the interac- 
tion network is self-organized. The parameter a determines the level number 
of hierarchical structures. For large a, a self-similar-like hierarchical structure 
is created, in which the population sizes and interaction ranges increase geo- 
metrically as the level is increased. The stationary solutions to Eq. (1) are not 
unique as shown in a globally coupled model [Hj, however, similar forms of hier- 
archical structures were obtained from different initial conditions. If the initial 
condition is completely localized as x.- L — 1 for i = N/2 and X4 = for the other 
i, the localized state is maintained in the time evolution of Eq. (1), since new 
populations cannot be born from zero. As a natural modified model, a diffusion 
process may be introduced in the population dynamics as 

= x i(Y*. w hj x j ~ c o) + D{x i+ i - 2xi + Xi-i), for i = 1, ■ • • ,N, (4) 
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where D is a diffusion constant. The diffusion process corresponds to the mi- 
gration in ecosystems. New populations are migrated at the sites where Xi — 
originally, and then the localized state becomes unstable. We have checked that 
similar hierarchical structures are naturally created from the localized initial 
conditions even for very small D. We have used a special form of interac- 
tions. However, the dynamical instability of the homogeneous state is due to 
the negative derivative of dwij(r)/dr at r — 1 and other forms of interactions 
with the same characteristics can induce similar hierarchical structures. The 
self-organization of the hierarchical structure owing to mutual interactions is 
suggestive of the formation of complex networks in ecosystems, economic and 
social systems [§]. 
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